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Heat capacity of matter is considered to be its most important property because it holds information about 
system's degrees of freedom as well as the regime in which the system operates, classical or quantum. Heat 
capacity is well understood in gases and solids but not in the third main state of matter, liquids, and is not 
discussed in physics textbooks as a result. The perceived difficulty is that interactions in a liquid are both 
strong and system- specific, implying that the energy strongly depends on the liquid type and that, therefore, 
liquid energy can not be calculated in general form. Here, we develop a phonon theory of liquids where this 
problem is avoided. The theory covers both classical and quantum regimes. We demonstrate good 
agreement of calculated and experimental heat capacity of 21 liquids, including noble, metallic, molecular 
and hydrogen-bonded network liquids in a wide range of temperature and pressure. 

Research into liquids has a long history starting from the same time when the theory of gases was developed, 
forming the basis for our current understanding of the gas state of matter 1 . Yet no theory of liquid heat 
capacity currently exists, contrary to gases or, for that matter, solids. Consequently, liquid heat capacity is 
not discussed in textbooks 1 " 6 , and presents a challenge at student lectures 7 . 

Liquids defy common approaches used to discuss the other two main states of matter, gases and solids. 
Interactions in a liquid are strong, therefore treating them as a small perturbation as is done in the theory of 
gases is not an option 2 . Atomic displacements in a liquid are large, therefore expanding the energy in terms of 
squared atomic displacements and considering the remaining terms small, as is done in the theory of solids, does 
not appear justified either. Strong interactions in a liquid appear to result in liquid energy being strongly 
dependent on the type of interactions, leading to a statement that liquid energy can not be calculated in general 
form 2 , in contrast to solids and gases where E = 3NT (k B = 1) and E= -NT at high temperature. 

Previous research into liquids presents an interesting case when long-lived absence of crucial experimental data 
resulted in using our human experience and intuition about liquids instead. Liquids have been viewed to occupy a 
state intermediate between gases and solids. Liquids flow, and share this fundamental property with gases rather 
than solids. As a result, liquids were mostly viewed as interacting gases, forming the basis for previous theoretical 
approaches i_ 6 . On the other hand, Frenkel noted that the density of liquids is only slightly different from that of 
solids, but is vastly different from gas densities i. Frenkel has subsequently made another important proposition 
which has further put liquids closer to solids in terms of their physical properties. He introduced liquid relaxation 
time, t, as the average time between two consecutive atomic jumps at one point in space, and immediately 
predicted the solid-like ability of liquids to support shear waves, with the only difference that a liquid does not 

support all shear waves as a solid does, but only those with frequency larger than - because the liquid structure 

T 

remains unaltered and solid-like during times shorter than t. Frenkel's ideas did not find due support at his time, 
partly because of the then existing dominating school of thought advocating different views 8 , and consequently 
were not theoretically developed. 

It took many years to verify Frenkel's prediction 9 " 13 . With the aid of powerful synchrotron radiation sources 
that became available fairly recently, and about 50 years after Frenkel's work, it has become apparent that even 
low- viscous liquids maintain high-frequency modes with wavelengths close to interatomic separations 13 . We note 
that being the eigenstates of topologically disordered liquids, these modes are phonon-like in a sense that they are 
not entirely harmonic vibrations as in crystals. 

The data from the experimental advances above, together with theoretical understanding of the phonon states 
in a liquid due to Frenkel, raise an important question of whether the phonon theory of liquid thermodynamics 
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Figure 1 | Experimental and calculated c^for liquid metals. Experimental c v are at ambient pressure, with electronic contribution subtracted 1819 . Values 
of t d used in the calculation are 0.95 ps (Cs), 0.27 ps (Ga), 0.49 ps (Hg), 0.31 ps (In), 0.57 ps (K), 0.42 ps (Na), 0.64 ps (Pb), 0.74 ps (Rb) and 0.18 ps (Sn). 
Values of Goo are 0.17 GPa(Cs),2.39 GPa(Ga),1.31 GPa(Hg),1.58 GPa(In),1.8 GPa(K),3.6 GPa(Na),1.42 GPa(Pb),0.25 GPa(Rb)and3 GPa(Sn). 
Experimental a (Ref. 19 ) are3-10" 4 K" 1 (Cs), 1.2-KT 4 K" 1 (Ga), 1.8-HT 4 K" 1 (Hg), 1.11* 10" 4 K" 1 (In), 2.9- 1(T 4 K" 1 (K), 2.57- 1(T 4 K" 1 (Na), 
3-10" 4 K" 1 (Pb),3'l(T 4 K- 1 (Rb) and 0.87* 10" 4 K" 1 (Sn). Values of a used in the calculation are 3.8* 10" 4 K" 1 (Cs), 1.2- 1(T 4 K" 1 (Ga), 1.6- 10" 4 K" 1 
(Hg), 1.25- 10" 4 K" 1 (In), 2.9-10 " 4 K" 1 (K), 2.57- 1(T 4 K" 1 (Na), 3*10" 4 K" 1 (Pb), 4.5- 10" 4 K" 1 (Rb) and 1.11-10" 4 K" 1 (Sn). 



could be constructed, similar to the phonon theory of solids. Below 
we develop a phonon theory of thermodynamics of liquids that cov- 
ers both classical and quantum regimes. We find good agreement 
between calculated and experimental heat capacity for many liquids 
in a wide temperature and pressure range. 

Results 

Theory. There are two types of atomic motion in a liquid: phonon 
motion that consists of one longitudinal mode and two transverse 

modes with frequency co > co F , where co F = — is Frenkel frequency, 

T 

and the diffusional motion due to an atom jumping between two 
equilibrium positions. In turn, the phonon and diffusional motion 
consists of kinetic and potential parts, giving the liquid energy as 

E = Ki + P l + K s (co>coi)+P s (co>co F ) + K d + P d (1) 

where Ki and Pi are kinetic and potential components of the 
longitudinal phonon energy, K s (w > w F ) and P s (w > w F ) are 
kinetic and potential components of the energy of shear phonons 
with frequency w > co F , and K d and P d are kinetic and potential 
energy of diffusing atoms. The energy of longitudinal mode is the 
same as in a solid, albeit different dissipation laws apply at low and 



high frequency 1 . In the next paragraph, we re- write E in the form 
convenient for subsequent calculations. 

P d can be omitted because P d <C P s (co > co F ). This is because the 
absence of shear modes with frequency w < co F due to the dynamic 
Frenkel mechanism implies smallness of low-frequency restoring 
forces and therefore smallness of low- frequency potential energy of 
shear modes: P s (w < co F ) <C P s (co > co F ), where P s (w < w F ) is the 
potential energy of low- frequency shear modes 1415 . Instead of low- 
frequency shear vibrations with potential energy P s (co < <x> F ), atoms 

2n 

in a liquid "slip" and undergo diffusing motions with frequency — 

T 

and associated potential energy P d . P d is due to the same interatomic 
forces that P s (co < ct) F ), giving P d ~ P s (co < co F ). Combining this with 
P s (w < co F ) <C P s (w > co F ), P d <^ P s (w > w F ) follows. Re-phrasing 
this, were P d large and comparable to P s (co > <x> F ), strong restoring 
forces at low frequency would result, and lead to the existence of low- 
frequency vibrations instead of diffusion. We also note that because 
Pi ~ P s , P d <^P s (cd> co F ) gives P d <C Pb further implying that P d can 
be omitted in Eq. (1). 

Then, E = K + P t + P s (w > cl> f ), where total kinetic energy in a 
3 

liquid, K (K = - N /c B T in the classical case), includes both vibrational 
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Figure 2 | Experimental and calculated c^for noble liquids. Experimental c v and f] are taken from the NIST database 17 at pressures 378 MPa (Ar), 
200 MPa (Kr), 70 MPa (Ne), and 700 MPa (Xe). Values of t d used in the calculation are 0.5 ps (Ar), 0.67 ps (Kr), 0.31 ps (Ne) and 0.76 ps (Xe). Values of 
Goo are 0.18 GPa (Ar), 0.23 GPa (Kr), 0.12 GPa (Ne) and 0.45 GPa (Xe). Experimental values of a calculated from the NIST database at the corresponding 
pressures above are 3.6- 1(T 4 K _1 (Kr), 7.7* 10~ 3 K" 1 (Ne) and 4.1 -10" 4 K" 1 (Xe). Values of a used in the calculation are 1.3* 10" 3 K" 1 (Kr), 6.8* 10" 3 K" 1 
(Ne) and 3.6* 10" 4 K -1 (Xe). For Ar, calculating c v in harmonic approximation (a = 0) gives good agreement with experimental c v . 



and diffusional components. E can be re-written using the virial 

. Ei , , , E s (w>(D?) 
theorem P/ = — and P s (w>w^) = (here, P and E refer 

to their average values) and by additionally noting that K is equal to 
kinetic energy of a solid and can therefore be written, using the virial 
theorem, as the sum of kinetic terms related to longitudinal and shear 

Ei E s E s E s ((d>cd?) „ 
waves: K = 1 , giving E = E\-\ 1 . Finally noting 

that E s = E s (co < co F ) + E s (co > co F ), where the two terms refer to 
their solid-state values, liquid energy becomes 



3*=*b+(l + ?)£ 



hcoi 



2 J ^ fh(Di\ 

1 expf — 1-1 



(4) 



t-i 17 I 17 / \ , E S(W<CJ ¥ ) 

E = Ei + E s {co><jo?) H 



(2) 



The energy of one longitudinal mode, the first term in Eq. (2), can 
be calculated by substituting the sum in Eq. (4), E, with Debye 
vibrational density of states for longitudinal phonons, 

g(co)= — a* 2 , where co D is Debye frequency. The normalization 

w D 

of g(w) reflects the fact that the number of longitudinal modes is 
N. Integrating from 0 to co D gives 2 = NTD(— ^j, where 



Each term in Eq. (2) can be calculated using the phonon free energy, 

F p h = E 0 + T In ^1 — exp ^ — ^ , where E 0 is the energy of 

zero -point vibrations 2 . In calculating the energy, 
dF t 



z 3 dz 



exp(z) — 1 



is Debye function 2 . The energy of two shear 



modes with frequency w > co F , the second term in Eq. (2), can be 
similarly calculated by substituting £ with density of states 
6N 



E ph =F ph -T^, we take into account the effect of thermal expan- sH=^co , where the normalization accounts for the number 



dT 



dwj 



sion, important in liquids. This implies — — # 0, contrary to the 

dT 

harmonic case, and gives 



(Dj — T 



d(Dj 

If 



(3) 



Using quasi-harmonic approximation Griineisen approximation 

dw occd T n 

gives — = , where a is the coefficient of thermal expansion . 

dT 2 

Putting it in Eq. (3) gives 



of shear modes of 2iV. Integrating from co F to co D gives 
Z = 2NId(*^>) _2Nr^ 3 D^. E s (co < o, F ) in the last 
term in Eq. (2) is obtained by integrating £ from 0 to w F with the 

same density of states, giving £ = 2NT \-—\ D ( ^— ^ J . Putting all 

\co D J \ T J 

terms in Eq. (4) and then Eq. (2) gives finally the liquid energy 

— K)H^)-£)X^)) <5 > 

where we have omitted zero-point vibration energy E 0 because below 
we calculate the derivative of E. We note that in general, E 0 is tem- 
perature-dependent because it is a function of co F and t (see Eq. (2)). 
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Figure 3 | Experimental and calculated c^for molecular liquids. Experimental c v and f] are taken from the NIST database 17 at pressures 50 MPa (CH 4 ), 
55 MPa (CO), 0.1 MPa (F 2 ), and 170 MPa (H 2 S), 65 MPa (N 2 ) and 45 MPa (0 2 ). Experimental c v has the rotational contribution subtracted (see text). 
Experimental rj for F 2 is from Ref. 26 . Values of t d used in the calculation are 0.44 ps (CH 4 ), 0.55 ps (CO), 0.71 ps (F 2 ), 0.65 ps (H 2 S), 0.38 ps (N 2 ) and 0.46 
ps (0 2 ). Values of Goo are 0.1 GPa (CH 4 ), 0.11 GPa (CO), 0.23 GPa (F 2 ), 0.15 GPa (H 2 S), 0.14 GPa (N 2 ) and 0.22 GPa (0 2 ). Experimental values of a 
calculated from the NIST database at the corresponding pressures above are 3.1 -10" 3 K" 1 (CH 4 ), 4.M0" 3 K" 1 (CO), 4.31- 10" 3 K" 1 (F 2 ), 1.2- 10" 3 K" 1 
(H 2 S), 3.9- 10" 3 K" 1 (N2) and 4.4- 10~ 3 K" 1 (0 2 ). Values of a used in the calculation are 1.4- 10 3 K _1 (CH 4 ), 5.5- 10" 3 K" 1 (CO), 4.31 • 10" 3 K" 1 (F 2 ), 
6.5- 10~ 4 K" 1 (H 2 S), 3.9- 10~ 3 K" 1 (N 2 ) and 4.4- 10~ 3 K" 1 (0 2 ). 



However, this becomes important at temperatures of several K only, 
whereas below we deal with significantly higher temperatures. 

Eq. (5) spans both classical and quantum regimes, an important 
feature for describing the behavior of liquids discussed below. 

We note that the presented phonon theory of liquids operates at 
the same level of approximation as Debye theory of solids, particu- 
larly relevant for disordered isotropic systems 2 such as glasses and 
liquids. The result for a harmonic solid follows from Eq. (5) when co F 
= 0, corresponding to infinite relaxation time, and a = 0. 

Comparison with experimental data. We now compare Eq. (5) to 
the experimental data of heat capacity. We have used the National 
Institute of Standards and Technology (NIST) database 17 that 
contains c v for many liquids, and have chosen monatomic nobel 
liquids, molecular liquids, as well as hydrogen-bonded network 
liquids. We have used c v measured on isobars. We aimed to check 
our theoretical predictions in a wide range of temperature, and 
therefore selected the data at pressures exceeding the critical 
pressures of the above systems where they exist in a liquid form in 
the broad temperature range. As a result, the temperature range in 



which we calculate c v is about 100-700 K for various liquids (see 
Figures 1, 2, 3 and 4). Experimental c v of metallic liquids were taken 
from Refs. 1819 at ambient pressure. We note that previously, we have 
considered c v of some of these metals on the basis classical approach 
only 14 . The total number of liquids considered is 21. 

We numerically evaluated Debye functions in Eq. (5) for each 

1 dE 

temperature and calculated c v = We have taken viscosity data 
from the NIST database at the same pressures as c v and converted it 

to t using the Maxwell relationship t = , where y\ is viscosity and 

2xc 2tcGoo 

Goo is infinite-frequency shear modulus, giving co F = — = • 

T f] 

Viscosities of liquid metals and F 2 were taken from Refs. 20-25 and 
Ref. 26 , respectively. 

We note that Debye model is not a good approximation in 
molecular and hydrogen-bonded systems where the frequency 
of intra-molecular vibrations considerably exceeds the rest of fre- 
quencies in the system (e.g. 2260 K in 0 2 , 3340 K in N 2 , 3572 K in 
CO). However, the intra-molecular modes are not excited in the 
temperature range of experimental c v (see Figures 3-4). Therefore, 
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Figure 4 | Experimental and calculated c v for H 2 0 and D 2 0. 

Experimental c v and rj are taken from the NIST database 17 at pressures 
150 MPa (H 2 0) and 80 MPa (D 2 0). Experimental c v has rotational and 
configurational contribution subtracted (see text). Values of t d used in the 
calculation are 0.15 ps (H 2 0) and 0.53 ps (D 2 0). For D 2 0, experimental a 
calculated from the NIST database at the corresponding pressure above is 
9.4- 10" 4 K _1 , and a used in the calculation is 1.1 -10" 3 K _1 . For H 2 0, 
calculating c v in harmonic approximation (a = 0) gives good agreement 
with to experimental c v . 



the contribution of intra-molecular motion to c v is purely rotational, 
c rot . The rotational motion is excited in the considered temperature 
range, and is therefore classical, giving c rot = R for linear molecules 



such as CO, F 2 , N 2 and 0 2 and c n 



3R 



for molecules with three 



rotation axes such as CH 4 , H 2 S, H 2 0 and D 2 0. Consequently, c v for 
molecular liquids shown in Figures 3-4 correspond to heat capacities 
per molecule, with c rot subtracted from the experimental data. In this 
case, N in Eq. (5) refers to the number of molecules. 

We also note that in H 2 0, approximately half of experimental c v is 
the configurational contribution due to the temperature-induced 
structural modification of the hydrogen-bonded network, although 
the precise value of this contribution is not currently settled 27,28 . The 
structural modification includes changing coordinations of H 2 0 
molecules during the continuous transition between the low- density 
and high-density liquid in the wide temperature range 28 . Similarly to 
the phonon theory of solids, effects of structural changes are not 
accounted for in the general theory presented here. Therefore, 

experimental c v for H 2 0 and D 2 0 shows c v with half subtracted from 

. 3R 

their values due to the configurational term and c = — subtracted 
further as discussed above. Due to the approximate way in which the 



configurational contribution is treated for H 2 0 and D 2 0, the agree- 
ment between the experimental and calculated c v should be viewed as 
qualitative, in contrast to the quantitative agreement for the rest of 
liquids considered. As for other molecular liquids above, the result- 
ing c v represents heat capacity per molecule. 

Experimental and calculated c v for 21 noble, metallic, molecular 
and hydrogen-bonded network liquids are shown in Figures 1-4. 
Figures 1-4 make one of the central points of our paper: the proposed 
phonon theory of liquids gives good agreement with experimental 
data. Importantly, calculated c v has no free adjustable parameters, 
but depends on <x> D , a and Goo (see Eq. (5)) which are fixed by system 

properties. <x> D , a and Goo that give the best agreement in Figures 1-4 

2tt 

are close to their typical values. t d = — used in the calculation (see 

the captions in Figures 1-4) are consistent with the known values 
that are typically in the 0.1-1 ps range 11,29 . For monatomic liquids, t d 
were taken as experimental t d in corresponding solids. Similarly, the 
difference between experimental a and a used in the calculation is 
within 30% on average. Finally, Goo used in the calculation are on the 
order of GPa typically measured 11,29 . 

Discussion 

We first note that some of the considered liquids are in the quantum 
regime at low temperature. Taking T as the lowest temperature in 

Figures 1-4 and co D = where t d are given in the captions of 

Figures 1-4, we find that — ^- is in the range 0.1-3 for various liquids. 

Consequently, some liquids can be described by classical approxi- 
mation fairly well, whereas others require quantum approach: for 

example, > 1 for Ne, 0 2 , N 2 , F 2 , CH 4 and CO, implying pro- 
gressive phonon excitation with temperature increase. 

We observe that c v decreases with temperature in Figures 1-4. 
There are two main competing effects that contribute to c v in Eq. 

(5). First, temperature increase results in the increase of cl> f = — • 

T 

Consequently, c v decreases as a result of the decreasing number of 
shear waves that contribute to liquid energy and c v . Second, c v 
increases with temperature due to progressive phonon excitation 
as discussed above. The first effect dominates in the considered tem- 
perature range, and the net effect is the decrease of c v with temper- 
ature seen in Figure 1. 

We further observe that c v changes from approximately 3 to 2, 
showing a universal trend across a wide range of liquids. This can be 

understood by noting that at low temperature, I — ) and the sec- 

\codJ 

ond term in the last bracket in Eq. (5) are small, giving c v close to 3 
(including the contribution from the anharmonic term). At high 
temperature when co F ~ co D and Debye functions are close to 1, 
the term in the second bracket in Eq. (5) is 2, giving c v ~ 2. 

We note that on further temperature increase, c v decreases from 2 
to its gas value of 1.5 17 . This is related to the progressive disappear- 
ance of longitudinal phonons, an effect discussed in our forthcoming 
paper. 

The decrease of c v seen in Figures 1-4 does not need to be generic. 
Indeed, when the phonon excitation dominates at low temperature, 
as in the case of strongly quantum liquids such as H 2 or He, c v 
increases with temperature 17 . 

Good agreement between theoretical and calculated c v makes sev- 
eral important points. First, it is interesting to revisit the statement 
that liquid energy can not be calculated in general form because 
interactions are both strong and system- specific 2 . In the proposed 
phonon theory of liquids, this difficulty does not arise because strong 
interactions are treated from the outset as in the theory of solids. 

Second, the proposed phonon theory of liquids has an advantage 
over the previous approach where liquid potential energy is 
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5 



calculated from correlation functions and interatomic interactions 2 ' 3 . 
Starting from the earlier proposals 30 , this approach was developed in 
several directions (see, e.g., Refs. 2,3 ' 31 " 34 ), and can be used to 
calculate the energy of simple liquids where interactions are pair 
and short-ranged and correlations are two-body as is the case in, 
for example, noble liquids or hard-sphere models. On the other hand, 
the calculations become intractable in general case 3,33 , particularly 
when interactions and correlations become complex (e.g., many- 
body) as in the liquids discussed above, precluding the calculation 
of c v and, consequently, understanding and interpreting experi- 
mental data. On the other hand, if the phonon states of the liquid 
depend on z only as proposed by Frenkel 1 , the liquid energy and c v 
depend implicitly on z only, even though correlation functions and 
interatomic interactions may affect both c v and z in a complex way. 
Then, the relationship between c v and z becomes fairly simple (see 
Eq. (5)) and explains the experimental behavior of a wide class of 
liquids, both simple and complex. 

The last assertion is important for a general outlook at liquids: 
despite their apparent complexity, understanding their ther- 
modynamics may be easier than previously thought. Indeed, we 
have good understanding of solid thermodynamics based on pho- 
nons no matter how complicated interactions or structural corre- 
lations in a solid are. Our current results suggest that the same can 
apply to liquids. 
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